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Abstract 

We present a simple and general method to generate a set of basis functions suitable for parametrizing 
the anisotropy of a given physical property in the presence of symmetry constraints. This mathematical 
tool has direct applications in various fields, for instance, to parametrize the orientation-dependence of 
interface energies or to represent the so-called constituent strain elastic energy of superlattice structures 
in the long wavelength limit. The proposed method can be easily implemented using generic linear 
algebra operations without having consider many different subcases that depend on the point group 
symmetry considered. The method exploits a direct correspondence between spherical harmonics and 
polynomial functions of a unit vectors expressed in tensor notation. The method has been implemented 
in the "genes" code of the Alloy Theoretic Automated Toolkit (ATAT). 

1 Introduction 

The expansion of a physical property in terms of a basis of known functions is a very common procedure in 
physics. In this paper, we focus on generating bases suitable for expanding a function on the unit sphere, 
that is, a function of direction describing the anisotropy of a given physical property. 

This mathematical tool has direct applications in various fields, for instance, to parametrize the orientation- 
dependence of interface energies [1], to represent the charge densities of rigid molecules in X-ray analysis 
[2] or to describe the orientation-dependence of the so-called constituent strain elastic energy of superlattice 
structures in the long wavelength limit [3,4]. Of course, it is well-known that, in the absence of symmetry 
constraints, spherical harmonics provide a complete solution to the problem. The problem becomes interest- 
ing and non-trivial when symmetry constraints reduce the number of degrees of freedom in the expansion, 
which is the topic of the present paper. One distinguishing feature of the proposed method is that it can be 
easily implemented using generic linear algebra operations without having consider many different subcases 
that depend on the point group considered. The algorithm proposed herein has been implemented in the 
"genes" code of the Alloy Theoretic Automated Toolkit (ATAT) [5-8]. The method also highlights a seldom 
mentioned characterization of spherical harmonics: They are simply polynomials in the components of a 
unit vector. 

2 Method Outline 

Polynomials are known to form a complete basis for any continuous function over a bounded region (e.g. 
the unit sphere). Hence, in particular, they form a complete basis for any continuous function defined over the 
surface of the unit sphere. Let u be a three-dimensional unit vector (e.g. u — (cos {9) sin (0) , sin {6) sin (0) , cos (0))). 
Any continuous function / of direction u can therefore be represented as 

oo 

f{u) = J2A(^^u^ (1) 

L=0 



1 



where we use the short-hand notation 




(with A^^^u^ defined as constant) and where A^^},,^^ is a rank L tensor that is symmetric under permutation 
of the indices (since permutations of the indices does not change the polynomial we can, without loss of 
generality, limit ourselves to such symmetric tensors). If the function / (u) is constrained by symmetry, such 
constraints can then be implemented by restricting the tensors A'-^^ to obey suitable invariance with respect 
to all symmetry operations in a given point group [9]. As explained in more detail in the Section 3, this 
can be simply accomplished by considering 3^ noncolinear trial tensors, and obtaining symmetrized tensors 
by averaging each trial tensor with all its transformations by each point group symmetry operation. The 
desired result is obtained after eliminating colinear symmetrized tensors. 

An additional step is needed because expansion (1) is a bit redundant, since the polynomial (w,f + + uf) 
is constant over the unit sphere. This would imply that nonzero coefficients A^^^ for L > could give rise to 
a constant / (u), which is undesirable. This can be avoided by projecting each A^^^ onto the space orthogonal 
to tensors giving rise to polynomials that can be factored as 

{ul+ul+ul)A^'^-^^u^-^ (2) 

for some tensor A^^~'^^ of rank L — 2. (It is not necessary to consider higher powers of (wf + w| + W3) in 
this factorization, because A^^~^^u^~^ could include additional (wf + + M3) factors as a special case.) A 
simple algorithm to accomplish the symmetrization and this projection is provided in Section 3. 

The result of this procedure is an expression for the tensor A'^^ as a sum of fc^ symmetry-constrained 
and non-redundant components A^^'''^: 

k=l 

where the tensors A'-^''^^ are fixed and determined by symmetry while the coefficients CL,k are completely 
unrestricted. Upon substitution into (1) we obtain a symmetry-constrained expansion: 

00 Kl 
L=Ofe=l 



3 Algorithm 

Let us first define a few convenient symbols. 



Let S denote a 3 x 3 matrix representing a point symmetry operation in Cartesian coordinates and let 
the corresponding function S (^A^^^) applied to a tensor A^^^ of rank L be defined as: 



1=1 iL = l k=l 

and let <S denote a set of such matrices that defines the point group of interest. 

• Let P denote a permutation vector (i.e. an L-dimensional vector containing all the number {1, . . . ,L} 
not necessarily in increasing order) and let the function P (•) be defined as: 



ip(i)f-tip(L) 



and let V^^^ denote the set of all such permutations for a given L. 



2 



• Let # denote the number of elements in a set. 

• Let C^^'^-', . . . , C'^^'^^^ be a set of rank L tensors defining linear constraints on the generated tensor 
basis, i.e. the A^^'''^ generated must be orthogonal to all C^^'"^\ according to the inner product 

3 3 
ii = l ii, = l 

(If M = 0, no constraints are imposed.) 

• Let vec (A^^)) denote a vectorization of the tensor A^^^ (i.e. a 3^-dimensional column vector containing 
all elements of the tensor A^^^) and let unvec (•) denote the reverse operation. 

Our algorithm for generating a basis for tensors of rank L obeying symmetric constraints (defined by 

a point group 5), indices permutation invariancc constraints (defined by the set V^^^) and some linear 
constraints (defined by a basis of tensors C*^^'^\ . . . , C*^^'^)) is then as follows: 

1. If M > 0, define B to be the nonzero and non linearly dependent columns of the matrix: 



B = 



(C(^'i)),...,vcc(c(^''^)) 



2. Set A; = 



3. Consider a set of distinct trial tensors ^(^'*) for t = 1, . . . , 3^, each consisting of a single element set 
to 1, with all remaining elements set to O."*^ 

(a) For each trial tensor A^^'*\ set 

Q = unvec {(j - B (B'^B) B^^ vec (a^^-*) 



(or simply Q = if M = 0). 

(b) Calculate 



(c) If Q is nonzero and (for A; > 0) not linearly dependent with the [A^^'^^ . . . jd^^'*^'] , then increment 
k and set = Q. 

4. Finally, set = k and orthogonalize (and normalize) the element of . . . A^^'^^^j using the 

Gram-Schmidt procedure. 

The harmonics of order up to L"^^ are then generated by calling the above routine for L = 1, . . . , L™*^, 
setting M = if i < 2 and otherwise setting the constraints C^^'''^ to be 

with C'ii,...,^^ = Si^ i^A^^~lf^ 

where the A^^~^'''^ are also generated with the above routine, called with rank L — 2, the same point group 
S, the permutation set and M = (no constraints C^^"^''')). 



^For added efficiency, one can limit the trial tensors to those whose nonzero element A.i^^,,,^i^ obeys ii < 12 < • • • < il- 
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4 Discussion 



4.1 Symmetrization technique 



It is instructive to see why the method used for symmetrization in the algorithm of Section 3 actually works. 
For an arbitrary trial tensor A^^'*^ we can verify that the symmetrized tensor 



Q 



ses 



obeys T (^Qj = Q for any operation T G S. Indeed, calculate 



ses 



where the second equality holds because T (S {■)) is just another operation in S {■) € S and two distinct 
S {■) cannot be mapped onto the same symmetry operation by applying T (•), since each element of a group 
admits an inverse. Since the symmetrization procedure is a linear projection, choosing the trial tensors A^^'*^ 
to be an orthogonal basis is sufficient to generate a basis for the space of symmetrized tensors. 

A similar argument holds for invariance under permutations P of the indices. Finally, note that applying 
a point group operation S {■) to a tensor A^^^ that is invariant to indices permutations yields a tensor with 
the same property: 



j{Pl),...,j{PL) 



3 3 L 

1 = 1 ii = l fe=l 



3 3 L 

ii = l iL = '^ k=l 

3 3 L 

H = l *i=l k=l 

3 3 L 

X/ ■ ■ ■ X/ ^n, n ^ikjk = 

ii = l ii,=l fe=l 



where we have used the fact re-ordering the sums or the product has no effect and the invariance of A^^^ under 
permutation. This shows that symmetrizing the tensor after making it invariant to indices permutations 
does not undo the permutation invariance. 



4.2 Relationship to spherical harmonics 

It is interesting to observe that expansion (3) coincides (apart from an inconsequential linear transformation) 
with spherical harmonics when no symmetry constraints arc imposed. The easiest way to see this is to 
compare (3) with the eigenfunction of the Schrodinger equation for some spherically symmetric potential 
selected so that the eigenfunctions involve polynomials. We can use any convenient radial potential because 
we only focus on the angular part. Consider a spherically symmetric harmonic potential, whose eigenstates 
are polynomials times a spherically symmetric Gaussian. The Gaussian is constant over the unit sphere, 
so we are left with only a polynomial as the angular dependence. Moreover, it is well-known that the 
order of that polynomial is equal to ni + n2 + n^, the sum of three principal quantum numbers of the 
harmonic oscillator along each dimension. This sum is also (up to a constant scaling and shift) the energy 
of the system. It follows that there is a direct correspondence between all terms in (3) sharing the same 
value of L and all eigenfunctions sharing the same energy. The different terms sharing the same L thus 
correspond to eigenstates with different angular momentum projections. We can verify that the number of 
terms (in the case of a spherically symmetric potential) matches the number of spherical harmonics for a 
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given value of L. Indeed, the number of distinct terms of total power L in a polynomial in v variables 



(^J^) ~ (2) = {L + 2){L + 1) /2—L (L — 1) /2 = 2L+1, exactly the number of spherical harmonics associated 
with angular momentum is L. Hence the dimension of the space spanned by the spherical harmonics for 
a given L is the same as the dimension of the space spanned by our polynomials. Both spaces include 
polynomials of order L on the unit sphere that are not colinear and it follows that both bases must span the 
same space. Hence both expansions, truncated to the same L, span the same space. 

4.3 Relationship to other approaches 

Suitable treatments for special cases already exist in the literature. Notable examples include the cubic 
harmonics (e.g. [10,11]) and harmonics for hexagonal symmetry (e.g. [12,13]). A very general treatment 
can be found in [2] . Although very complete, this treatment does not lends itself to a simple implementation: 
the point group and its orientation has to be identified, not just as a list of symmetry operations, but 
recognized by name as one of the known point groups, so that one can lookup the specific rules applying 
to that point group. Based on the point group category found (e.g., cubic or hexagonal), a superset of 
harmonics is selected. Then, based on the specific point group found, "index rules" are applied to eliminate 
those harmonics that should vanish by symmetry. This treatment is ideally suited for researchers wanting 
to manually construct a basis based on the knowledge of the point group, as the different case are nicely 
classified by point group. However, a computer program implementing the method would also necessarily 
contain a large number of tests and subcases. It would also have to rotate the symmetries into a standard 
"setting" to use the tabulated index rules. Moreover, if one wishes to handle other points group that are 
not special cases of cubic or hexagonal symmetries (e.g. icosahedral symmetry or other noncrystallographic 
point groups), a different superset of harmonics and index rules must be constructed. 

In contrast, the approach proposed herein works directly with the symmetry operation in matrix form 
(which are easy to determine) and no classification into categories of point group is needed. The possibility 
of having point group in different orientation (or "settings" ) is automatically handled, with no extra coding 
effort. The algorithm only relies on basic linear algebra operations and handles any point groups, not just 
those for which supersets of harmonics have already been constructed. The proposed method is related to 
the one proposed in [14] to generate tensor bases, although additional steps, provided herein, were needed to 
formally show that such tensors bases can be used to generate direction-dependent harmonics and to avoid 
redundant harmonics via a projection scheme. 

5 Examples 

Figure 1 shows the harmonics generated by the proposed algorithm for each of the crystallographic point 
groups. The coefficients A^-'^^ in electronic form can be found in the Supplement Material (in the ATAT 
format described in Section 6, along with the input files needed to generate the harmonics). Figure 2 shows 
the result of a similar exercise for selected noncrystallographic point groups. 

6 Using the genes code 

The code takes, as an input, either a point group (specified via generators) or a structural information from 
which it determines the point symmetry automatically. It outputs the harmonics in a file, in form that is easy 
to read into a computer code and outputs the harmonics in human-readable form on the standard output. 
Some of the file formats contain extraneous items not needed for harmonic generation per se (marked in 
italics below), but that are included to ensure compatibility with other portions of the ATAT package. 

6.1 Input files 

By default, the code reads in structural information (from the lat . in file by default — a alternate file name 
can be specified with the -1 option) and determines the point group automatically. The lat. in file has the 
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Figure 1: Symmetry-adapted spherical harmonics for each crystahographic point groups (up to L=6). Group 
1 (with 48 harmonics) is omitted. The inset shows the location of the cartesian axes on the unit sphere (where 
directions of highest symmetry are alligned, whenever possible) and the color scheme used to represent the 
function. Unique axis (when appropriate) is chosen to be z. 
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Figure 2: Symmetry-adapted spherical harmonics for selected noncrystallographic point groups (up to L=6), 
represented using the same conventions as in Figure 1. 



following format. 

1. First, the coordinate system a,b,c is specified, either as 

[a] [b] [c] [a] m [7] 

or in terms of Cartesian coordinates, one axis per line: 

[qx] [ay] [a^] 
[b.] [by] [b.] 

[Cx] [Cy] [C,] 

2. Then the lattice vectors u, v, w are listed, one per line, expressed in the coordinate system just defined: 

[Ua] [Ub] [Uc] 

[Va] [Vb] [Vc] 

[Wa] [Wb] [Wc] 

3. Finally, the position Xi and type(s) tii,t2i, ■ ■ ■ of atom for site i are given, expressed in the same 
coordinate system as the lattice vectors: 

[Xal] [Xbl] [Xcl] [tll,t21---] 
[Xa2] [Xb2] [Xcl] [il2,i22,---] 



An example of such file, for an Al-Ti alloy adopting the hep crystal structure is: 



3, 


.1 3.1 5.062 90 90 120 


(Coordinate system: a h c. q B 7 notation) 


1 





(Primitive unit cell: one vector per line 





1 


expressed in multiples of the above coordinate 





1 


system vectors) 





Al.Ti 


(Atoms in the lattice) 


0. 


.6666666 0.3333333 0.5 Al.Ti 





As an alternative to providing the above structural information, the user can provide generators of the 
point group (which could also be the whole point group) and the code will complete the full point group 
automatically. This input file is called sym.in and has the format: 

[Number symmetry operations given as generators] 

[3x3 matrix representing a symmetry operation in Cartesian coordinates] 


[another 3x3 matrix representing a symmetry operation in Cartesian coordinates] 




6.2 Output file 

The harmonics are output in the file harm. out, which has the following format: 
For each harmonic: 

1 




tensor 
[rank] 

[3's repeated 'rank' times] 

rank J 

[the 3''™''^ elements of the tensor] 

In addition, the standard output displays the harmonics in human-readable format, with x, y, z denoting 
the components of a unit vector. 
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6.3 Command line options 

-l=[string] Input file defining tlie lattice (default: lat.in) 

-s Specify point group directly via generator (in sym.in) 

-z=[real] Tolerance for finding symmetry operations (default: 2e-4) 

-r=:[int] Maximum Rank of the harmonic 

-mr=[int] Minimum Rank of the harmonic (default 1) 

-sig=[int] Number of significant digits printed (default: 5) 

Invoking the code without any options displays help. At the minimum, the user must specify the -r 
option. All other options are optional. 

Acknowledgements 

This work is supported by XSEDE computing resources and by the National Science Foundation under Grant 
No. DMR-0907669. 

References 

[1] M. Siegel, M. J. Miksis, and P. W. Voorhees. Evolution of material voids for liighly anisotropic surface 
energy. J. Mech. Phys. Sol, 52:1319—1353, 2004. 

[2] M. Kara and K. Kurki-Suonio. Symmetrized multipole analysis of orientational distributions. Acta 
Cryst., A37:201, 1981. 

[3] D. B. Laks, L. G. Fcrreira, S. Froyen, and A. Zunger. Efficient cluster expansion for substitutional 

systems. Physical Review B, 46:12587, 1992. 

[4] V. Ozolins, C. Wolverton, and Alex Zunger. Effects of anharmonic strain on the phase stability of 
epitaxial films and super lattices: Applications to noble metals. Phys. Rev. B, 57:4816, 1998. 

[5] A. van de Walle and G. Ceder. Automating first-principles pliase diagram calculations. Journal of 
Phase Equilibria, 23:348, 2002. 

[6] A. van de Walle, M. Asta, and G. Ceder. The Alloy Theoretic Automated Toolkit: A user guide. 
CALPHAD Journal, 26:539, 2002. 

[7] A. van de Walle. Multicomponent multisublattice alloys, nonconfigurational entropy and other additions 
to the Alloy Theoretic Automated Toolkit. Calphad Journal, 33:266, 2009. 

[8] A. van de Walle. The alloy theoretic automated toolkit, 2001. http://www.its.caltech.edu/~avdw/atat/. 

[9] J. F. Nye. Physical Properties of Crystals: Their Representation by Tensors and Matrices. Oxford 

University Press, USA, 1985. 

[10] S. L. Altmann and A. P. Cracknell. Lattice harmonics I. cubic groups. Rev. Mod. Phys., 37:19, 1965. 

[11] J. Muggli. Cubic harmonics as linear combinations of spherical harmonics. J. of Applied Mathematics 
and Physics, 23, 1972. 

[12] S. L. Altmann and C. J. Bradley. Lattice harmonics II. hexagonal close-packed groups. Rev. Mod. 
Phys., 37:33, 1965. 

[13] C. Varney and G. L. W. Hart. A coherency strain model for hexagonal-close-packed alloys. TMS Letters, 
1:35, 2004. 

[14] A. van de Walle. A complete representation of structure-property relationships in crystals. Nature 
Materials, 7:455, 2008. 



8 



